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Using a modified Lennard- Jones model for anisotropic particles, we present results of molecu- 
lar dynamics simulation in two dimensions. In one-component systems, we find crystallization, a 
Berezinskii-Kosterlitz-Thouless phase, and a structural phase transition, as the temperatures is low- 
ered. In the lowest temperature range, the crystal is composed of three martensitic variants on 
a hexagonal lattice, exhibiting the shape memory effect. With addition of larger spherical parti- 
cles (impurities), these domains are finely divided, yielding glass with slow time evolution. With 
increasing the impurity size, the structural or translational disorder is also proliferated. 

PACS numbers: 81.30.Kf, 61.43.Fs, 61.72.-y, 64.70. kj 



Certain anisotropic particles such as KCN form a cu- 
bic crystal and, at lower temperatures, they undergo an 
order-disorder phase transition, where the crystal struc- 
ture changes to a noncubic one. Furthermore, with addi- 
tion of impurities, the phase ordering often occurs only 
on small spatial scales, where heterogeneous orientation 
fluctuations are pinned [T]. In such systems softening 
of the shear modulus is observed, indicating direct cou- 
pling between the molecular orientation and the acoustic 
phonons, and the molecules often have dipolar moments, 
yielding dielectric anomaly. These systems with frozen 
disorder have been identified as orientational glass. As 
a similar example, metallic ferroelectric glass, called re- 
laxor, with frozen polar nanodomains have been studied 
extensively [2 . Recently, a system of off-stoichiometric 
intermetallic Ti-Ni was shown to be glassy martensite or 
strain glass, exhibiting the shape-memory effect and the 
superelasticity 5] ■ For a- one-component system of hard 
spheroids, Frenkel and Mulder ^4^ performed Monte Carlo 
simulation to find isotropic liquid, nematic liquid, orien- 
tationally ordered solid, and orientationally disordered 
(plastic) solid. Theoretical approaches on strain glass so 
far have been a phase-field theory with elastic field and 
a random temperature [S] and a spin-glass theory with 
elastic long-range interaction [6j. 

In this Letter, we propose a simple microscopic model 
exhibiting orientation-martensitic phase transitions and 
glass behavior. In two dimensions, we suppose ellip- 
tic particles interacting via an angle-dependent Lennard- 
Jones potential, where their positions are and their ori- 
entation vectors are rii = (cos 6';, sin 6*.;) [i = I,-- - ,A^). 
There can be two particle species 1 and 2 with radii 
CTi and (72 and numbers A'^i and N2. We set N = 
Ni + N2 — 4096 and change the composition c — N2/N. 
The pair potential Uij between particles i € a and j S /3 
(a, /3 = 1, 2) is expressed as 



e is the characteristic interaction energy. The particle 
anisotropy is taken into account by the angle factor. 



(2) 



where f, 



^rij represents the direction of r.y. We 



introduce the anisotropy parameters xi and X2 for the 
two species. We truncate the above potential for rij > Vc 
with the cut-off being Vc — Sci. We also set Cij = Uij 
at rij = Vc to ensure the continuity of the potential at 
'''ij — Tci so there remains a weak angle-dependence in 
Cij. Our potential is analogous to the Gay-Berne poten- 
tial for anisotropic molecules [7], which has been used to 
simulate mesophases of liquid crystals, and the Shintani- 
Tanaka potential with five-fold symmetry yielding frus- 
trated particle configurations [8]. 

The total potential and kinetic energies are U = 
Y.^<JU^, and K = J2^[m\dr,/dt\^ + Ia\d9^/dt\^]/2, re- 
spectively, where the two species have a common mass 
m and inertia momenta Ii and l2- The Newton equa- 
tions of motions are 
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where r 



Since we treat equilibrium or nearly steady states at a 
given temperature T, we attach a Nose-Hoover thermo- 
stat [9 to all the particles by adding the thermostat terms 
in Eq.(3). Space, time, and T will be measured in units 
of (Ti, To = (Ti-^/m/e, and e/fcs, respectively. In our sim- 
ulation, we started with a liquid at T = 2, quenched 
the system to T = 0.35 below the melting temperature 
Tjn ~ 1.0, and annealed it for 9000to. We then lowered 
T to a final low temperature. 

Assuming that the particles of the second species are 
spherical and larger, we set X2 = 0, 02lo\ = 1.2 or 1.4, 
and c = 0,0.1 or 0.2. From Eq.(2) the particles of the 
species 1 have short and long diameters given by — 
2^l^a\ and = (1 -|- Vlxij^^^as^ so their molecular area 

- Cij , 

(1) is 5*1 ~ Tiasai/A and their inertia momentum is Ii — 
{a} + al)m/4:, while I2 = cr|m/2 and S2 = ■n2^^/^al. 
\i'ij\^ <^ai3 = (cq + o',3)/2, and The packing fraction (A^i^i -I- N2S2)/V is fixed at 0.95 
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FIG. 1. Orientation angle 9i (left) and order parameter am- 
plitude (right) for c = and xi = 0.1 at T = 0.074, 
0.07, and 0.05 from above. Bottom left: Structure factor 
SQ{k) of the orientation fluctuations, growing for small k in 
the range T2 < T < Ti. Bottom right: Average amplitude 
(g^) = E, li/N vs r for xi = 0.1, 0.2, 0.3, and 0.4. 



and the system length is about 70(Ti. For each particle 
i of the first species, we introduce the orientation tensor 

Qi = {QifiA (Mi ^ ^x,y) as 

Q, = (1 + iV^)-i(n,:n, + ^ n-jUj)- 7 /2 

J G bonded 



g,(d,d,- / /2), 



(4) 
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where / is the unit tensor and di is the director with 
\di\ = 1. The summation is over the bonded particles 



FIG. 2. Shape memory effect under uniaxial stretching along 
the y axis at T = 0.02 for c = and Xi = 0.1. Left: Strain 
£ VS applied stress {oyy) in units of For {(Jyy) > 0.075, 

there remains only the variant elongated along the y axis. 
After this cycle, the residual strain vanishes upon heating to 
T — 0.1. Right: Fractions of the three variants during the 
cycle, which are stretched along the three crystal axes. 



(|rjj| < SiJi) of the first species with being the num- 
ber of these bonded particles. When a hexagonal lattice 
is formed, it includes the second nearest neighbor parti- 
cles. The angle of di varies more smoothly than 9i. The 
amplitude q, is given by qf = 2 J2f,,u Qtfuy- 

First, we show numerical results in the one-component 
case (c = 0) with xi =0.1 to study the orientation phase 
transition on a hexagonal lattice. Here we use the peri- 
odic boundary condition at fixed volume, but essentially 
the same results followed at zero pressure. In Fig.l, we 
show the orientation angle 9i of all the particles (left) 
and the order parameter amplitude qf (right) at T = 
0.074,0.07, and 0.05. From the angle snapshots we rec- 
ognize emergence of three variants with lowering T due to 
the underlying hexagonal lattice. The left bottom panel 
shows the structure factor SQ{k) = (IQjfcl^) ^"-"^ ~ 
J2jiQjxx — Q jyy) Gyip{ik ■ r j) , while the right bottom 
panel displays the average (9^) = over all the 

particles for xi = 0.1,0.2,0.3, and 0.4. The orientation 
order develops gradually in a narrow region T2 < T < Ti, 
where T2 - 0.070 and Ti - 0.076 for xi = 0.1. The 
Ti and T2 increase with increasing xi- In this temper- 
ature window, a Berezinskii-Kosterlitz-Thouless (BKT) 
phase [TOl [11] is realized between the low-temperature 
martensitic phase and the high-temperature orientation- 
ally disordered phase, where the orientation fluctuations 
are much enhanced at long wavelengths. Though our sys- 
tem size is still small, SQ{k) apparently grows as k^~'^ 
for k ^ 0.5, where rj depends on T (where 77 = 0.05 at 
T = 0.074). We should note that Bates and Frenkel [12] 
performed Monte Carlo simulation of two-dimensional 
rods to find the Kosterlitz-Thouless phase transition. 

For T < T2, the three variants become distinct with 
sharp interfaces. The surface tension between the vari- 
ants is about O.le/af for xi = 0.1 (and is about 0.2e/cr^ 
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for xi = 0.2). In the pattern at T = 0.05 in Fig.l, the 
junction angles, at which two or more domain bound- 
aries intersect, are muUiples of 7r/6. This geometrical 
constraint stops the domain growth at a characteristic 
size even without impurities '13J. Similar patterns were 
observed on hexagonal planes in a number of experiments 
[2] and were reproduced by phase-field simulation \TE\ . 

In our model, the orientation order induces lattice de- 
formations. As a result, softening of the shear modulus fi 
occurs near the transition [T], while the bulk modulus K 
remains of order 20e/af. In fact, for c = and xi — 0.1, 
we have /x ~ 3 at T = 0.1, ^ 1.5 at T = 0.08, and 
/i ^ 5 for T < 0.05 in units of e/crf. Each variant at 
low T is composed of isosceles triangles elongated along 
one of the crystal axes, where side lengths are 1.21(71 and 
l.llcTi for xi = 0.1 at T = 0.05 in Fig.l. 

In our system, there arises a shape memory effect [3]. 
In Fig. 2, we applied a stress {<7yy) along the y axis at 
T = 0.02 [16], treating the surface along it as a free 
boundary. Initially, the fractions of the three variants 
were nearly close to 1/3 and one variant was elongated 
along the y axis. For very small £ < 2 x 10"'^, the 
system deformed elastically with /z ~ 2. However, for 
2 X lO^'^ < £ < 0.075, the fraction of the favored vari- 
ant increased up to unity, while those of the disfavored 
ones decreased. This inter-variant transformation occurs 
without defect formation. In the next step (ayy) was de- 
creased slowly from O.le/ai. On this return path, the 
solid was composed of the favored variant only. At van- 
ishing stress, there remained a remnant strain, but it 
disappeared upon heating to T = 0.1 above the transi- 
tion. Here, we may define the effective shear modulus 
^le by [de/d{ayy)]-^ = ^K^lJ{K + ^^e) = 4^e. Then 
/ie ^ 0.1 — 0.8 during the inter- variant transformation 
and /Xg ^ 5 on the return path. 

Next, in Fig. 3, we present examples of strain glass with 
impurities, where <JiI<J\ = 1.2 and xi = 0.1. For c = 0.1 
and 0.2, there appeared a few tens of particles with co- 
ordination numbers different from six in a single crystal. 
In our model, the elliptic particles tend to be parallel to 
the surface of the larger spherical ones, resulting in an- 
choring of the orientation. We can see that the size of the 
domains decreases with increasing c, where the impuri- 
ties suppress the development of the orientation order. In 
addition, the BKT phase disappeared in these examples. 
We also observed a shape-memory effect even in orien- 
tational glass states |3j, where small disfavored domains 
were replaced by the favored ones upon stretching. 

Now, we discuss the dynamics. Let us consider the 
time-dependent angle-distribution function, 

G{tM = + to) - e.ito) - ^))/iVi, (5) 

j 

where the average (• • •) is taken over the initial time to 
and over several runs. We are interested in the first two 
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FIG. 3. Frozen patterns of angle 6i (left) and order parameter 
amplitude qi (right) with impurities for c = 0.1 and 0.2, where 
Xi ~ 0.2 at T = 0.05. Bottom left: Expanded snapshot of 
Qidi in a box in the upper panel, showing pinned mesoscopic 
order or strain. Bottom right: (g^) vs T for various c. 
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FIG. 4. Orientation relaxation time ri from the time- 
correlation function Gi{t) for (a) c = and xi = 0.1, (b) 
c = 0.1 and xi = 0.1, (c) c = and Xi = 0.2, and (d) c = 0.1 
and xi ~ 0.2. It represents the turnover time. For (a) and (c), 
Ti grows steeply in the Berezinskii-Kosterlitz-Thouless phase 
{T2 < T < Ti). For strain glass (b) and (d), this phase is 
nonexistent and ri grows as T is lowered. 
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c=0.20, Xi=0.2, a2/ai = 1.4 




FIG. 5. Orientation angle 9i (left) and six-fold bond orienta- 
tion angle Oi in Eq.(7) (right) in polycrystal for a2/o"i = 1-4 
and c = 0.2, where Xi = 0.2 and T = 0.05. 



moments Gi(t) and G2{t). For p = 1, 2, • • • , we define 

<.27r 

Gp(t)= / d^G{t, if) cos{p^)/2n, (6) 
Jo 

which decays from unity on a time scale of Tp. Here, ri 
is the inverse frequency of the turnover motions 9j {t) — >■ 
9j (t) ± TT at low T, while T2 is the randomization time of 
cos(29j{t) — 26j{0)). Each turnover motion takes place 
quickly. We find that G{t, ip) exhibits a peak of the form 
A{t) exp[—{ip — 7r)^/2(T^]/-\/27rcr for small (p — n at low T, 
where a ^ 0.45. The coefficient A{t) grows linearly as 
O.St/ri for t Ti and tends to a constant (= 1/2) for 
i > Ti. The fitting Gi(i) = exp[-(i/Ti)^] fairly holds, 
where /3 decreases from unity to about 0.5 as T is lowered. 

For c = 0, Ti increases steeply in the BKT phase and 
Tj"^ cx cxp{-To/T) for T < T2, where (a) Tq - 0.40 
at xi = 0.1 and (c) Tq - 1.2 at xi = 0.2 in Fig.4. In 
addition, n - T2 for T > Ti but ts/ti > 1 for T < Tj. In 
fact, for xi = 0.1, the ratio is about 10^ at T = 0.07 and 
is about 10^ at T = 0.06. On the other hand, in glassy 
states with impurities, the relaxation behavior is more 
complicated due to the pinning effect, but the turnover 
motions still occur and T2 ^ ri holds. Figure 4 shows 
that Ti for c = 0.1 is longer in the disordered phase but 
is shorter in the ordered phase than in the pure system. 

For (T2/o'i ~ 1-2 and for c = 0.1 and 0.2, the crystal 
structure is little affected by the orientation fluctuations. 
For a larger size ratio, the structural or positional disor- 
der is more enhanced, eventually resulting in polycrystal 
and glass [17]. In Fig. 5, we reafize a polycrystal state 
for (T2/cri = 1.4, xi = 0.2, and c = 0.2, where black 
points represent the impurities. The left panel displays 
9j , where there remains noticeable orientation order with 
(g^) = 0.2. The right panel displays the positional six- 
fold orientation angle aj [TTJ[T7]. Here, for each elliptic 



particle j, we define aj in the range < < 7r/3 by 

exp[6i9jk] = Zj exp[6iaj], (7) 

fctEbondcd 

where 9jk is the angle of rjk = Vk — rj with respect to 
the X axis. We set \rjk\ < 1.7iTi and Zj > 0. 

In summary, we have presented an angle-dependent 
Lennard- Jones potential to simulate orientation or 
martensitic transitions. We have added impurities, which 
pin orientation and strain fluctuations on mesoscopic 
scales. In future, we should examine the impurity pin- 
ning on the glass transition in detail by systematically 
changing the composition and the size ratio [17) . Com- 
petition of the orientational and translational glass be- 
haviors should also be studied. We will shortly report 
three-dimensional simulation results, where inclusion of 
the dipolar interaction will enrich the problem. 
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